Code
import geopandas as gpd
#PASDA libraries shapefile
libraries = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PhiladelphiaLibraries201302/PhiladelphiaLibraries201302.shp" )
#Swiming Pools
pools = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Swimming_Pools/PPR_Swimming_Pools.shp" )
#Playgrounds and Rec Centers
playgrounds = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Playgrounds/PPR_Playgrounds.shp" )
#Parks
parks = gpd.read_file("/Users/macbook/Desktop/Python/Week 14/PPR_Program_Sites.geojson" )
Code
for name, gdf in {
"libraries" : libraries,
"pools" : pools,
"playgrounds" : playgrounds,
"parks" : parks
}.items():
print (name, gdf.crs)
libraries EPSG:2272
pools EPSG:3857
playgrounds EPSG:3857
parks EPSG:4326
Code
TARGET_CRS = "EPSG:2272"
layers = {
"libraries" : libraries,
"pools" : pools,
"playgrounds" : playgrounds,
"parks" : parks
}
for name, gdf in layers.items():
if gdf.crs != TARGET_CRS:
layers[name] = gdf.to_crs(TARGET_CRS)
libraries = layers["libraries" ]
pools = layers["pools" ]
playgrounds = layers["playgrounds" ]
parks = layers["parks" ]
Code
libraries["type" ] = "Library"
pools["type" ] = "Pool"
playgrounds["type" ] = "Playground"
parks["type" ] = "Park"
Code
import pandas as pd
import geopandas as gpd
import folium
import censusdata
Code
u18_vars = [
"B01001_001E" , # total population
"B01001_003E" ,"B01001_004E" ,"B01001_005E" ,"B01001_006E" ,
"B01001_027E" ,"B01001_028E" ,"B01001_029E" ,"B01001_030E" ,
]
acs_u18 = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
u18_vars
).reset_index()
Code
edu_vars = [
"B15003_001E" , # total 25+
"B15003_022E" ,"B15003_023E" ,"B15003_024E" ,"B15003_025E"
]
acs_edu = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
edu_vars
).reset_index()
Code
inc_vars = ["B19013_001E" ]
acs_inc = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
inc_vars
).reset_index()
Code
def make_geoid(df):
df["GEOID" ] = df["index" ].apply (
lambda x: "" .join([part[1 ] for part in x.geo])
)
return df
Code
acs_u18 = make_geoid(acs_u18)
acs_edu = make_geoid(acs_edu)
acs_inc = make_geoid(acs_inc)
Code
acs_u18["under_18" ] = acs_u18[
[
"B01001_003E" ,"B01001_004E" ,"B01001_005E" ,"B01001_006E" ,
"B01001_027E" ,"B01001_028E" ,"B01001_029E" ,"B01001_030E"
]
].sum (axis= 1 )
acs_u18["pct_under_18" ] = (
acs_u18["under_18" ] / acs_u18["B01001_001E" ] * 100
)
Code
acs_edu["bachelors_plus" ] = acs_edu[
["B15003_022E" ,"B15003_023E" ,"B15003_024E" ,"B15003_025E" ]
].sum (axis= 1 )
acs_edu["pct_bachelors_plus" ] = (
acs_edu["bachelors_plus" ] / acs_edu["B15003_001E" ] * 100
)
Code
acs_inc = acs_inc.rename(columns= {
"B19013_001E" : "median_income"
})
Code
acs_full = (
acs_u18[["GEOID" , "pct_under_18" ]]
.merge(
acs_edu[["GEOID" , "pct_bachelors_plus" ]],
on= "GEOID" ,
how= "left"
)
.merge(
acs_inc[["GEOID" , "median_income" ]],
on= "GEOID" ,
how= "left"
)
)
Code
import geopandas as gpd
# Read tracts geometry and merge the consolidated ACS table (acs_full)
tracts_geom = gpd.read_file(
"https://www2.census.gov/geo/tiger/TIGER2022/TRACT/tl_2022_42_tract.zip"
)
tracts_geom = tracts_geom[tracts_geom["COUNTYFP" ] == "101" ]
tracts = tracts_geom.merge(
acs_full,
on= "GEOID" ,
how= "left" ,
)
Code
# --- social infrastructure attribute renames and combine ---
libraries = libraries.rename(columns= {"ASSET_NAME" : "name" })
pools = pools.rename(columns= {"amenity_na" : "name" })
playgrounds = playgrounds.rename(columns= {"park_name" : "name" })
parks = parks.rename(columns= {"park_name" : "name" })
Code
social_infra = gpd.GeoDataFrame(
pd.concat([
libraries[["name" , "type" , "geometry" ]],
pools[["name" , "type" , "geometry" ]],
playgrounds[["name" , "type" , "geometry" ]],
parks[["name" , "type" , "geometry" ]],
], ignore_index= True ),
crs= libraries.crs
)
Code
social_infra_web = social_infra.to_crs("EPSG:4326" )
Code
# Project tracts and points to a projected CRS
tracts_proj = tracts.to_crs("EPSG:2272" )
social_infra_proj = social_infra.to_crs("EPSG:2272" )
Code
infra_join = gpd.sjoin(
social_infra_proj,
tracts_proj[["GEOID" , "geometry" ]],
how= "left" ,
predicate= "within"
)
Code
infra_counts = (
infra_join
.groupby("GEOID" )
.size()
.reset_index(name= "facility_count" )
)
Code
infra_counts.head()
infra_counts["facility_count" ].describe()
count 225.000000
mean 3.382222
std 2.032250
min 1.000000
25% 2.000000
50% 3.000000
75% 4.000000
max 13.000000
Name: facility_count, dtype: float64
Code
type_cols = {
"Library" : "library_count" ,
"Playground" : "rec_center_count" ,
"Rec Center" : "rec_center_count" ,
"Park" : "park_count" ,
"Pool" : "pool_count"
}
infra_by_type = (
infra_join
.groupby(["GEOID" , "type" ])
.size()
.unstack(fill_value= 0 )
.reset_index()
.rename(columns= type_cols)
)
for col in type_cols.values():
if col not in infra_by_type.columns:
infra_by_type[col] = 0
Code
tracts_access = (
tracts_proj
.merge(infra_counts, on= "GEOID" , how= "left" )
.merge(infra_by_type, on= "GEOID" , how= "left" )
)
# Replace NA with 0
count_cols = ["facility_count" , "library_count" , "rec_center_count" , "park_count" , "pool_count" ]
tracts_access[count_cols] = tracts_access[count_cols].fillna(0 )
Code
tracts_access_web = tracts_access.to_crs("EPSG:4326" )
Code
count_field_options = [
("Social Infrastructure Count:" , ["facility_count" ]),
("Libraries:" , ["library_count" , "Library" ]),
("Rec Centers:" , ["rec_center_count" , "Rec Center" , "Playground" ]),
("Parks:" , ["park_count" , "Park" ]),
("Pools:" , ["pool_count" , "Pool" ]),
]
tooltip_fields = []
tooltip_aliases = []
for alias, options in count_field_options:
for col in options:
if col in tracts_access_web.columns:
tooltip_fields.append(col)
tooltip_aliases.append(alias)
break
if not tooltip_fields:
tooltip_fields = ["facility_count" ]
tooltip_aliases = ["Social Infrastructure Count:" ]
Code
from folium.plugins import DualMap
dual_map = DualMap(
location= [39.9526 , - 75.1652 ],
zoom_start= 11 ,
tiles= "CartoDB positron"
)
Code
import numpy as np
def style_facilities(feature):
v = feature["properties" ].get("facility_count" )
if v is None or v == 0 or (isinstance (v, float ) and np.isnan(v)):
return {
"fillColor" : "#f0f0f0" ,
"color" : "#cccccc" ,
"weight" : 0.3 ,
"fillOpacity" : 0.45
}
return {
"fillColor" : cmap_facilities(v),
"color" : "#ffffff" ,
"weight" : 0.3 ,
"fillOpacity" : 0.9
}
Code
tracts_access_web["facility_count" ].describe()
count 408.000000
mean 1.865196
std 2.260431
min 0.000000
25% 0.000000
50% 1.000000
75% 3.000000
max 13.000000
Name: facility_count, dtype: float64
Code
import branca.colormap as cm
bins = tracts_access_web["facility_count" ].quantile(
[0 , 0.2 , 0.4 , 0.6 , 0.8 , 1 ]
).values
cmap_facilities = cm.StepColormap(
colors= cm.linear.Blues_06.colors,
index= bins,
vmin= bins[0 ],
vmax= bins[- 1 ]
)
Code
import censusdata
import geopandas as gpd
import pandas as pd
import folium
Code
acs_vars = {
"total_pop" : "B01001_001E" ,
"male_under_18" : [
"B01001_003E" ,"B01001_004E" ,"B01001_005E" ,"B01001_006E"
],
"female_under_18" : [
"B01001_027E" ,"B01001_028E" ,"B01001_029E" ,"B01001_030E"
]
}
tracts_acs = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
[acs_vars["total_pop" ]] +
acs_vars["male_under_18" ] +
acs_vars["female_under_18" ]
)
Code
tracts_acs["under_18" ] = tracts_acs[
[
"B01001_003E" ,
"B01001_004E" ,
"B01001_005E" ,
"B01001_006E" ,
"B01001_027E" ,
"B01001_028E" ,
"B01001_029E" ,
"B01001_030E" ,
]
].sum (axis= 1 )
Code
edu_vars = [
"B15003_001E" , # total 25+
"B15003_022E" , # bachelor's
"B15003_023E" , # master's
"B15003_024E" , # professional
"B15003_025E" # doctorate
]
Code
income_vars = ["B19013_001E" ]
Code
tracts_acs_edu = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
edu_vars
)
tracts_acs_inc = censusdata.download(
"acs5" ,
2022 ,
censusdata.censusgeo([("state" ,"42" ), ("county" ,"101" ), ("tract" ,"*" )]),
income_vars
)
Code
tracts_acs_edu = tracts_acs_edu.reset_index()
tracts_acs_inc = tracts_acs_inc.reset_index()
Code
tracts_acs_edu = tracts_acs_edu.rename(columns= {
"B15003_001E" : "total_25plus" ,
"B15003_022E" : "bachelors" ,
"B15003_023E" : "masters" ,
"B15003_024E" : "professional" ,
"B15003_025E" : "doctorate"
})
Code
# If total population is still B01001_001E, rename it first
tracts_acs = tracts_acs.rename(columns= {
"B01001_001E" : "total_pop"
})
# Recreate under_18 if needed
if "under_18" not in tracts_acs.columns:
tracts_acs["under_18" ] = tracts_acs[
[
"B01001_003E" ,"B01001_004E" ,"B01001_005E" ,"B01001_006E" ,
"B01001_027E" ,"B01001_028E" ,"B01001_029E" ,"B01001_030E"
]
].sum (axis= 1 )
# NOW create the percentage
tracts_acs["pct_under_18" ] = (
tracts_acs["under_18" ] / tracts_acs["total_pop" ] * 100
)
Code
if "bachelors" in tracts_acs_edu.columns and "pct_bachelors" not in tracts_acs_edu.columns:
tracts_acs_edu["pct_bachelors" ] = (
tracts_acs_edu["bachelors" ]
/ tracts_acs_edu["total_25plus" ] * 100
)
Code
tracts_acs_inc = tracts_acs_inc.rename(columns= {
"B19013_001E" : "median_income"
})
Code
tracts_acs_edu["GEOID" ] = tracts_acs_edu["index" ].apply (
lambda x: "" .join([part[1 ] for part in x.geo])
)
Code
tracts_acs_inc["GEOID" ] = tracts_acs_inc["index" ].apply (
lambda x: "" .join([part[1 ] for part in x.geo])
)
Code
tracts_acs["GEOID" ] = tracts_acs.index.map (
lambda x: (
x.geo[0 ][1 ] + # state
x.geo[1 ][1 ] + # county
x.geo[2 ][1 ] # tract
)
)
Code
tracts_acs[["total_pop" , "under_18" , "pct_under_18" ]].head()
Census Tract 1.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000101
1947
40
2.054443
Census Tract 1.02; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000102
2897
162
5.591992
Census Tract 2; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000200
3486
414
11.876076
Census Tract 3; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000300
3914
280
7.153807
Census Tract 4.01; Philadelphia County; Pennsylvania: Summary level: 140, state:42> county:101> tract:000401
2675
127
4.747664
Code
acs_full = (
tracts_acs[["GEOID" , "pct_under_18" ]]
.merge(
tracts_acs_edu[["GEOID" , "pct_bachelors" ]],
on= "GEOID" ,
how= "left"
)
.merge(
tracts_acs_inc[["GEOID" , "median_income" ]],
on= "GEOID" ,
how= "left"
)
)
Code
tracts = tracts_geom.merge(
acs_full,
on= "GEOID" ,
how= "left"
)
tracts_web = tracts.to_crs("EPSG:4326" )
Code
import folium
m_acs = folium.Map(
location= [39.9526 , - 75.1652 ],
zoom_start= 11 ,
tiles= "CartoDB positron"
)
Code
import numpy as np
bins_income = tracts_web["median_income" ].quantile(
[0 , 0.2 , 0.4 , 0.6 , 0.8 , 1 ]
).values
Code
import branca.colormap as cm
cmap_income = cm.StepColormap(
colors= cm.linear.YlOrRd_06.colors,
vmin= bins_income[0 ],
vmax= bins_income[- 1 ],
index= bins_income
)
Code
bins_bach = tracts_access_web["pct_bachelors_plus" ].quantile(
[0 , .2 , .4 , .6 , .8 , 1 ]
).values
Code
import branca.colormap as cm
import numpy as np
cmap_bach = cm.StepColormap(
colors= cm.linear.Purples_06.colors,
index= bins_bach,
vmin= bins_bach[0 ],
vmax= bins_bach[- 1 ]
)
cmap_bach.caption = "% Bachelor's Degree or Higher"
Code
def style_bachelors(feature):
v = feature["properties" ].get("pct_bachelors_plus" )
if v is None or np.isnan(v):
return {
"fillColor" : "#A9A9A9" ,
"color" : "white" ,
"weight" : 0.4 ,
"fillOpacity" : 0.85
}
return {
"fillColor" : cmap_bach(v),
"color" : "white" ,
"weight" : 0.4 ,
"fillOpacity" : 0.85
}
Code
folium.GeoJson(
tracts_web,
name= "% Bachelor’s Degree or Higher" ,
style_function= style_bachelors
).add_to(m_acs)
<folium.features.GeoJson at 0x14ee4e5d0>
Code
cmap_income = cm.linear.YlOrRd_09.scale(
tracts_web["median_income" ].min (),
tracts_web["median_income" ].max ()
)
Code
import numpy as np
bins_income = tracts_web["median_income" ].quantile(
[0 , 0.2 , 0.4 , 0.6 , 0.8 , 1 ]
).values
Code
import branca.colormap as cm
cmap_income = cm.StepColormap(
colors= cm.linear.YlOrRd_06.colors,
vmin= bins_income[0 ],
vmax= bins_income[- 1 ],
index= bins_income
)
Code
def build_acs_layer(name, field, palette):
values = tracts_access_web[field]
bins = values.quantile([0 , 0.2 , 0.4 , 0.6 , 0.8 , 1 ]).values
cmap = cm.StepColormap(
colors= palette,
index= bins,
vmin= bins[0 ],
vmax= bins[- 1 ]
)
cmap.caption = name
def style_fn(feature, field= field, cmap= cmap):
val = feature["properties" ].get(field)
if val is None or (isinstance (val, float ) and np.isnan(val)):
return {"fillColor" : "#d9d9d9" , "color" : "white" , "weight" : 0.3 , "fillOpacity" : 0.4 }
return {"fillColor" : cmap(val), "color" : "white" , "weight" : 0.3 , "fillOpacity" : 0.85 }
return {"name" : name, "field" : field, "style" : style_fn, "cmap" : cmap}
acs_layers = [
build_acs_layer("% Bachelor's Degree or Higher" , "pct_bachelors_plus" , cm.linear.Purples_06.colors),
build_acs_layer("% Under Age 18" , "pct_under_18" , cm.linear.Greens_06.colors),
build_acs_layer("Median Household Income" , "median_income" , cm.linear.YlOrRd_06.colors),
]
FINAL MAP
Code
import folium
from folium.plugins import DualMap
from folium.features import GeoJsonTooltip
import branca.colormap as cm
import numpy as np
Code
#Dual Map Setup
from folium.plugins import DualMap
import folium
from folium.features import GeoJsonTooltip
dual_map_edu = DualMap(
location= [39.9526 , - 75.1652 ],
zoom_start= 11 ,
tiles= "CartoDB positron"
)
Code
folium.GeoJson(
tracts_access_web,
style_function= style_facilities,
tooltip= GeoJsonTooltip(
fields= tooltip_fields,
aliases= tooltip_aliases,
localize= True ,
sticky= True
)
).add_to(dual_map_edu.m1)
# Legend (optional but fine here)
cmap_facilities.add_to(dual_map_edu.m1)
0.0 2.2 4.3 6.5 8.7 10.8 13.0
Code
for layer in acs_layers:
fg = folium.FeatureGroup(name= layer["name" ], show= (layer["field" ] == "pct_bachelors_plus" ))
folium.GeoJson(
tracts_access_web,
style_function= layer["style" ],
tooltip= GeoJsonTooltip(
fields= [layer["field" ]],
aliases= [f" { layer['name' ]} :" ],
localize= True ,
sticky= True
)
).add_to(fg)
fg.add_to(dual_map_edu.m2)
layer["cmap" ].add_to(dual_map_edu.m2)
folium.LayerControl(collapsed= False ).add_to(dual_map_edu.m2)
0.5299417064122947 16.5 32.5 48.4 64.4 80.3 96.29629629629629 % Bachelor's Degree or Higher
Code
0
42
101
026301
42101026301
263.01
Census Tract 263.01
G5020
S
406340
0
...
-075.1637161
POLYGON ((-75.16899 40.07147, -75.16864 40.071...
22.435737
28.359827
57992
0.0
0.0
0.0
0.0
0.0
1
42
101
029200
42101029200
292
Census Tract 292
G5020
S
1745920
28444
...
-075.1017831
POLYGON ((-75.11288 40.02746, -75.1128 40.0274...
23.318995
19.735431
82784
0.0
0.0
0.0
0.0
0.0
2
42
101
024400
42101024400
244
Census Tract 244
G5020
S
421189
0
...
-075.1638925
POLYGON ((-75.16912 40.02387, -75.16843 40.024...
33.115265
40.127077
56912
3.0
0.0
1.0
2.0
0.0
3
42
101
033200
42101033200
332
Census Tract 332
G5020
S
842716
0
...
-075.0449684
POLYGON ((-75.05464 40.04465, -75.05396 40.045...
29.119272
26.290400
80409
2.0
0.0
1.0
0.0
1.0
4
42
101
980200
42101980200
9802
Census Tract 9802
G5020
S
4916641
250846
...
-075.0443487
POLYGON ((-75.07448 40.08685, -75.07443 40.087...
6.666667
24.867725
-666666666
5.0
0.0
1.0
4.0
0.0
5 rows × 21 columns
Code
Make this Notebook Trusted to load map: File -> Trust Notebook
Code
folium.GeoJson(
tracts_access_web,
style_function= style_facilities,
tooltip= GeoJsonTooltip(
fields= tooltip_fields,
aliases= tooltip_aliases,
localize= True ,
sticky= True
)
).add_to(dual_map.m1)
<folium.features.GeoJson at 0x16ba98150>
Code
for layer in acs_layers:
fg = folium.FeatureGroup(name= layer["name" ], show= (layer["field" ] == "pct_bachelors_plus" ))
folium.GeoJson(
tracts_access_web,
style_function= layer["style" ],
tooltip= GeoJsonTooltip(
fields= [layer["field" ]],
aliases= [f" { layer['name' ]} :" ],
localize= True ,
sticky= True
)
).add_to(fg)
fg.add_to(dual_map.m2)
layer["cmap" ].add_to(dual_map.m2)
folium.LayerControl(collapsed= False ).add_to(dual_map.m2)
<folium.features.GeoJson at 0x16b7a7d10>
Code
# Add aligned legends for both map panes
cmap_facilities.add_to(dual_map.m1)
0.5299417064122947 16.5 32.5 48.4 64.4 80.3 96.29629629629629 % Bachelor's Degree or Higher
Code
Make this Notebook Trusted to load map: File -> Trust Notebook